Primero cargamos los paquetes que vamos a utilizar

library(tidyverse)
library(knitr)
library(here)

También los paquetes específicos para EDA, pero esta vez no lo mostramos en el documento final porque en el chunk ponemos echo = FALSE:

Y alguno más para modelos:

library(broom)
library(ggfortify)
library(GGally)
library(modelsummary)   #- remotes::install_github('vincentarelbundock/modelsummary')
library(equatiomatic)   #- remotes::install_github("datalorax/equatiomatic")
library(cluster)        
library(factoextra)   
library(FactoMineR)
library(rpart)
library(rattle)

0. Intro

Hemos visto en tutoriales anteriores como cargar, arreglar y manejar datos con R à la tidyverse. En este tutorial vamos a ver algunos paquetes y funciones que nos pueden ayudar a acelerar el análisis exploratorio inicial de nuestros datos.

R es un lenguaje/entorno que fue creado para hacer estadística y casi no hemos visto nada de modelos estadísticos. Al final de este tutorial veremos como estimar algunos modelos con R. Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo (solo algo) de Machine Learning.

Una de las máximas de un data scientist es, como señalan aquí, “Know Your Data”. Además nos explican que:

Data exploration is the art of looking at your data, rapidly generating hypotheses, quickly testing them, then repeating again and again and again. The goal of data exploration is to generate many promising leads that you can later explore in more depth.

With exploratory data analysis, you’ll combine visualisation and transformation with your curiosity and scepticism to ask and answer interesting questions about data

Simplemente dos recursos:

1. Datos

1.1 Datos originales

Son datos de bebes valencianos para los años 2010 y 2017. Los datos provienen del INE.Concretamente están aquí.

Cargamos los datos y el diccionario en memoria de R (Global Environment)

df_orig <- rio::import(here::here("datos", "nacidos", "df_partos_23v_CV_10_17.rds")) #- 90.236 x 23
diccionario <- rio::import(here::here("datos", "nacidos", "dicc_partos_23v_CV_10_17.xlsx"))

Los datos originales tienen 23 variables de 90236 partos.

1.2 Selección de muestra

Vamos a restringir el análisis a partos de un sólo bebe

df <- df_orig %>% filter(MULTSEX %in% c(1,2))  #- solo partos con 1 bebe (88.176)

Con lo que ahora tenemos 88176 registros

Además voy a quitar algunas variables:

df <- df %>% select(MUNPAR.n, SEMANAS, CESAREA.f, PESON1, SEXO1.f, NUMHV, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f, PAISNACM.n, PAISNACP.n, fecha_parto, fecha_parto_anterior)   #- me quedo con 14 variables

Finalmente, nuestro df tiene 14 variables de 88176 partos.


2. Exploración de los datos (EDA)

En este apartado haremos repaso de algunas funciones y paquetes que nos pueden ser de utilidad en el análisis exploratorio inicial.

2.1 Nombres de las variables

Conviene saber acceder a los nombres de las variables. Por ejemplo con base::names()

names(df)
#>  [1] "MUNPAR.n"             "SEMANAS"              "CESAREA.f"           
#>  [4] "PESON1"               "SEXO1.f"              "NUMHV"               
#>  [7] "EDADM.2"              "EDADP.2"              "ESTUDIOM.f"          
#> [10] "ESTUDIOP.f"           "PAISNACM.n"           "PAISNACP.n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

Y saber cambiarlo si lo necesitamos

names(df)[2] <- "nuevo_nombre_variable_2"
names(df)
#>  [1] "MUNPAR.n"                "nuevo_nombre_variable_2"
#>  [3] "CESAREA.f"               "PESON1"                 
#>  [5] "SEXO1.f"                 "NUMHV"                  
#>  [7] "EDADM.2"                 "EDADP.2"                
#>  [9] "ESTUDIOM.f"              "ESTUDIOP.f"             
#> [11] "PAISNACM.n"              "PAISNACP.n"             
#> [13] "fecha_parto"             "fecha_parto_anterior"

Ejercicio: por favor, volved a dejar el nombre de la segunda variable como estaba: SEMANAS

names(df)[2] <- "SEMANAS"
names(df)
#>  [1] "MUNPAR.n"             "SEMANAS"              "CESAREA.f"           
#>  [4] "PESON1"               "SEXO1.f"              "NUMHV"               
#>  [7] "EDADM.2"              "EDADP.2"              "ESTUDIOM.f"          
#> [10] "ESTUDIOP.f"           "PAISNACM.n"           "PAISNACP.n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

Si necesitamos arreglar los nombres, la función janitor::clean_names() es fantástica para esto. Lo hace automáticamente, además tiene opciones para hacerlo como más te guste.

df_jugar <- df
janitor::clean_names(df_jugar) %>% names()
#>  [1] "munpar_n"             "semanas"              "cesarea_f"           
#>  [4] "peson1"               "sexo1_f"              "numhv"               
#>  [7] "edadm_2"              "edadp_2"              "estudiom_f"          
#> [10] "estudiop_f"           "paisnacm_n"           "paisnacp_n"          
#> [13] "fecha_parto"          "fecha_parto_anterior"

2.2 Estructura y tipo de variables

Siempre-siempre hay que saber de que tipo son nuestras variables. Esto se puede hacer de muchas formas. Os muestro algunas:

str(df)            #- muestra la estructura interna de un objeto R
#> Classes 'tbl_df', 'tbl' and 'data.frame':    88176 obs. of  14 variables:
#>  $ MUNPAR.n            : chr  "Villajoyosa/Vila Joiosa, la" "Alzira" "Gandia" "Valencia" ...
#>  $ SEMANAS             : num  39 40 41 34 40 38 40 40 38 37 ...
#>  $ CESAREA.f           : Factor w/ 2 levels "Sin cesárea",..: 1 2 1 2 1 1 1 1 1 1 ...
#>  $ PESON1              : num  3000 3260 3650 1720 2800 ...
#>  $ SEXO1.f             : Factor w/ 2 levels "Varón","Mujer": 1 1 1 2 1 1 1 2 2 1 ...
#>  $ NUMHV               : num  0 0 0 2 2 0 0 0 0 0 ...
#>  $ EDADM.2             : num  34.4 21.3 24.6 35.2 40.9 ...
#>  $ EDADP.2             : num  NA NA NA NA NA NA NA NA NA NA ...
#>  $ ESTUDIOM.f          : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#>  $ ESTUDIOP.f          : Factor w/ 3 levels "Primarios","Medios",..: NA NA NA NA NA NA NA NA NA NA ...
#>  $ PAISNACM.n          : chr  "Bulgaria" "Bulgaria" "Bulgaria" "Bulgaria" ...
#>  $ PAISNACP.n          : chr  NA NA NA NA ...
#>  $ fecha_parto         : POSIXct, format: "2017-05-01" "2017-07-01" ...
#>  $ fecha_parto_anterior: POSIXct, format: NA NA ...
inspectdf::inspect_types(df)   #- muestra de que tipo son las variables
#> # A tibble: 4 x 4
#>   type             cnt  pcnt col_name 
#>   <chr>          <int> <dbl> <list>   
#> 1 numeric            5  35.7 <chr [5]>
#> 2 factor             4  28.6 <chr [4]>
#> 3 character          3  21.4 <chr [3]>
#> 4 POSIXct POSIXt     2  14.3 <chr [2]>

2.3 Estadísticos básicos de las variables

Yo suelo utilizar estas dos funciones

df_aa <- pjp.funs::my_df_estadisticos_basicos(df)  #- estadisticos basicos del df
df_bb <- pjp.funs::my_df_valores_unicos(df)        #- valores únicos de cada variable de df

La función summarytools::dfSummary() da un resumen muy útil

#library(summarytools)
view(dfSummary(df))     #- genera un fichero con un resumen muy útil y agradable de ver para cada variable

Otra alternativa es usar skimr::skim(df)

skimr::skim(df)

2.4 NA’s

Siempre hay que ver la existencia de NA’s en nuestro df. Tres paquetes que nos ayudan en esta tarea son DataExplorer, visdat y naniar.

DataExplorer::plot_missing(df)
naniar::gg_miss_var(df, show_pct = TRUE)                     #- nº de % de NA's en cada columna

El paquete visdat es útil para ver los NA’s pero el siguiente chunk no nos va a funcionar. A ver si conseguís que funcione:

visdat::vis_dat(df)

Como resulta que visdat se queja de que hay muchas variables, entonces, vamos a visualizar sólo las variables que sí tienen NA’s:

Ejercicio: seleccionar las variables que tienen algún NA.

zz_con_NAs  <- df %>% select_if(anyNA)
zz_factores <- df %>% select_if(is.factor)
zz_x        <- df %>% select(PESON1, EDADM.2, EDADP.2, SEXO1.f)

Visualicemos las variables con NA's

visdat::vis_dat(zz_con_NAs)

  • Otra forma de visualizar la co-ocurrencia de NA’s en las variables:
naniar::gg_miss_upset(df)

  • Más de NA’s con naniar: Concretamente veremos la ocurrencia de NA’s para los grupos/categorías de la variable ESTUDIOM.f1
naniar::gg_miss_var(df, facet = ESTUDIOM.f, show_pct = TRUE) #- faceted por la variable ESTUDIOM.f

  • Vamos a borrar df’s que no vamos a utilizar
rm(list = ls(pattern = "^zz"))  #- borra los objetos cuyo nombre empiece por zz
# rm(zz_con_NAs, zz_factores, zz_x)

¿Que hacemos con los NA’s?

Pues no está claro, depende … pero vamos a dejar un df SIN NAs; y vamos a aprovechar para practicar algo de manejo de datos con R:

Ejercicio: quita los NA’s del df

zz_1 <- df[complete.cases(df), ]             #- con base-R 
zz_2 <- df %>% filter(complete.cases(.))     #- con dplyr y base-R
zz_3 <- df %>% tidyr::drop_na()              #- con tidyr

Ok, si quitamos todos los registros/partos que tienen algún NA en alguna de las 14 variables nos quedamos con 26492 registros.

PERO, la última variable, fecha_parto_anterior tiene muchos NA’s porque para muchas mujeres ese parto/bebe es el primero, así que vamos a refinar la eliminación de NA’s

zz_4 <- df %>% tidyr::drop_na(-fecha_parto_anterior)    #- quito filas con NA sin tener en cuenta fecha_parto_anterior
zz_5 <- df %>% tidyr::drop_na(1:13 )                    #- quito filas con NA en las 13 primeras variables

Sí, vamos a hacerlo:

df <- df %>% tidyr::drop_na(-fecha_parto_anterior) 
rm(list = ls(pattern = "^zz"))  #- borra los objetos cuyo nombre empiece por zz

Ahora nuestro df tiene 57051 partos, registros o filas.

2.5 Prevalencia de un valor

A veces es útil esta función

visdat::vis_expect(df, ~.x == 1)  #- ver la prevalencia de un determinado valor (puede ser de utilidad, mira la ayuda!!)

valores_a_buscar <- c("Con cesárea", "Universidad")
visdat::vis_expect(df, ~.x %in% valores_a_buscar) 

Ejercicio: seleccionar los partos en los que ambos padres tienen estudios universitarios y el parto se realizó con cesárea.

zz <- df %>% filter(ESTUDIOM.f == "Universidad" & 
                    ESTUDIOP.f == "Universidad" &
                    CESAREA.f == "Con cesárea" )


2.6 variables numéricas

Vamos a analizar un poco las variables numéricas.

Ejercicio: crea un df solo con las variables numéricas

df_numeric <- df %>% select_if(is.numeric)

Por ejemplo, summarytools::descr() nos ayuda a calcular rápidamente estadísticos descriptivos de las variables numéricas.

summarytools::descr(df, style = "rmarkdown")

Descriptive Statistics

df

N: 57051

  EDADM.2 EDADP.2 NUMHV PESON1 SEMANAS
Mean 32.30 34.98 0.62 3237.90 38.96
Std.Dev 5.27 5.80 0.80 508.86 1.89
Min 13.01 14.59 0.00 450.00 20.00
Q1 29.19 31.44 0.00 2960.00 38.00
Median 32.69 34.86 0.00 3250.00 39.00
Q3 35.94 38.36 1.00 3560.00 40.00
Max 56.04 71.46 11.00 6350.00 45.00
MAD 4.95 5.19 0.00 444.78 1.48
IQR 6.75 6.93 1.00 600.00 2.00
CV 0.16 0.17 1.29 0.16 0.05
Skewness -0.38 0.25 2.06 -0.61 -2.31
SE.Skewness 0.01 0.01 0.01 0.01 0.01
Kurtosis 0.15 1.11 9.60 2.35 10.78
N.Valid 57051.00 57051.00 57051.00 57051.00 57051.00
Pct.Valid 100.00 100.00 100.00 100.00 100.00
zz <- summarytools::descr(df)


Y DataExplorer lo vamos a utilizar para obtener histogramas o gráficos de densidad para las variables numéricas

DataExplorer::plot_histogram(df)

DataExplorer::plot_density(df)


Matriz de correlaciones

Veamos la matriz de correlaciones entre las variables numéricas:

df %>% select_if(is.numeric) %>% corrr::correlate() %>% gt::gt()
SEMANAS PESON1 NUMHV EDADM.2 EDADP.2
SEMANAS NA 0.52092275 -0.05170781 -0.04590809 -0.02587726
PESON1 0.52092275 NA 0.06702148 0.00018377 0.02165019
NUMHV -0.05170781 0.06702148 NA 0.23533870 0.23102964
EDADM.2 -0.04590809 0.00018377 0.23533870 NA 0.64879388
EDADP.2 -0.02587726 0.02165019 0.23102964 0.64879388 NA

Otra forma de calcular y visualizar las correlaciones entre las variables de un df

DataExplorer::plot_correlation(df)

  • solo para las variables continuas
DataExplorer::plot_correlation(df, type = 'c') #- solo las continuas

Otra forma de visualizar las correlaciones:

zz <- df %>% inspectdf::inspect_cor() 
df %>% inspectdf::inspect_cor() %>% inspectdf::show_plot()

Boxplot frente a una v. categórica

Para ver diferencias en la distribución de las variables numéricas frente a una v. categórica es muy útil DataExplorer::plot_boxplot(). Por ejemplo frente a la variable ESTUDIOM.f:

DataExplorer::plot_boxplot(df, by = "ESTUDIOM.f")

O, por ejemplo frente a la variable CESAREA.f

my_vv <- names(df)[3]
my_vv <- "CESAREA.f"
DataExplorer::plot_boxplot(df, by = my_vv)

  • boxplot para las variables numéricas frente a 2 categóricas
df %>% explore::explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)

df %>% select(EDADM.2, PESON1,  ESTUDIOM.f, CESAREA.f) %>% explore::explore_all(target = CESAREA.f)


Scatterplots entre variables numéricas

DataExplorer::plot_scatterplot() nos permite hacer rápidamente un scatterplot de todas las variables del df frente a una variable, por ejemplo el peso del bebe: PESON1

my_vv <- "PESON1"
DataExplorer::plot_scatterplot(df, by = my_vv, sampled_rows = 1000L)

2.7 Variables categóricas

Si necesitásemos un df con todas las variables categóricas, ¿cómo lo obtendrías?

df_categoricas <- df %>% select_if(negate(is.numeric))     #- !!!!

Para visualizar rápidamente los valores de las variables categóricas, tenemos inspectdf::inspect_cat(), que nos devuelve un gráfico con la distribución de todas las variables categóricas.

inspectdf::inspect_cat(df) %>% inspectdf::show_plot(high_cardinality = 1)


  • o con DataExplorer::plot_bar()
DataExplorer::plot_bar(df)




2.8 PRUEBAS (pkgs con Shiny)

  • El paquete GGally: una extensión de ggplot2 que sirve para muchas cosas, por ejemplo:
df_jugar <- df %>% select(EDADM.2, SEXO1.f, SEMANAS, PESON1)
GGally::ggpairs(df_jugar, progress = FALSE)
GGally::ggscatmat(df_jugar, columns = c("EDADM.2", "SEMANAS", "PESON1"))

burro attempts to make EDA accessible to a larger audience by exposing datasets as a simple Shiny App

#- devtools::install_github("laderast/burro")
burro::explore_data(df, outcome_var = colnames(df))


Instead of learning a lot of R syntax before you can explore data, the explore package enables you to have instant success. You can start with just one function - explore() - and learn other R syntax later step by step

  • Podemos crear un shiny para explorar df:
explore(df)
  • o explorar una o varias variables.
df %>% explore(CESAREA.f)
df %>% explore(EDADM.2, ESTUDIOM.f, target = CESAREA.f)




Funciones útiles para limpieza

  • calcular las ocurrencias de un determinado valor
zz <- df %>% janitor::get_dupes(PESON1)    #- janitor::get_dupes() calcula los valores de PESON1 repetidos
  • quitar filas y columnas vacías
zz <- df %>% janitor::remove_empty(c("rows", "cols"))    #- quita variables y filas vacias
  • La función dplyr::distinct() es muy útil para ver las filas con valores únicos. Por ejemplo:
zz <- df %>% distinct(ESTUDIOM.f, CESAREA.f)
zz %>% kable()
ESTUDIOM.f CESAREA.f
Primarios Con cesárea
Universidad Sin cesárea
Primarios Sin cesárea
Medios Sin cesárea
Universidad Con cesárea
Medios Con cesárea
zz %>% kable(format = "html") %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))
ESTUDIOM.f CESAREA.f
Primarios Con cesárea
Universidad Sin cesárea
Primarios Sin cesárea
Medios Sin cesárea
Universidad Con cesárea
Medios Con cesárea
  • Este es el chunk favorito de David Robinson:
df %>% count(PAISNACM.n) %>% 
  mutate(PAISNACM.n = fct_reorder(PAISNACM.n, n)) %>% 
  arrange(desc(n)) %>% slice(1:10) %>% 
     ggplot(aes(PAISNACM.n, n)) + geom_col() + coord_flip() +
       labs(title = "Nº de bebes por nacionalidad de la madre")



3. Tablas

Muchas veces hay que presentar resultados básicos como una tabla de casos o de frecuencias. para ello tenemos muchas posibilidades:

3.1 Tablas con summarytools

El paquete summarytools.

summarytools is an R package providing tools to neatly and quickly summarize data. It can also make R a little easier to learn and to use, especially for data cleaning and preliminary analysis.

  • La función freq() provee tablas con conteos y frecuencias
summarytools::freq(df$SEXO1.f, style = "rmarkdown")

Frequencies

df$SEXO1.f

Type: Factor

  Freq % Valid % Valid Cum. % Total % Total Cum.
Varón 29534 51.77 51.77 51.77 51.77
Mujer 27517 48.23 100.00 48.23 100.00
<NA> 0 0.00 100.00
Total 57051 100.00 100.00 100.00 100.00
  • tabulación cruzada entre dos variables categóricas:
summarytools::ctable(df$SEXO1.f, df$CESAREA.f)

Cross-Tabulation, Row Proportions
SEXO1.f * CESAREA.f
Data Frame: df

CESAREA.f Sin cesárea Con cesárea Total
SEXO1.f
Varón 20472 (69.3%) 9062 (30.7%) 29534 (100.0%)
Mujer 19486 (70.8%) 8031 (29.2%) 27517 (100.0%)
Total 39958 (70.0%) 17093 (30.0%) 57051 (100.0%)

Incluso se puede hacer un test chi-cuadrado

summarytools::ctable(df$ESTUDIOM.f, df$CESAREA.f, chisq = TRUE)

Cross-Tabulation, Row Proportions
ESTUDIOM.f * CESAREA.f
Data Frame: df

CESAREA.f Sin cesárea Con cesárea Total
ESTUDIOM.f
Primarios 15734 (73.2%) 5769 (26.8%) 21503 (100.0%)
Medios 11764 (68.8%) 5335 (31.2%) 17099 (100.0%)
Universidad 12460 (67.5%) 5989 (32.5%) 18449 (100.0%)
Total 39958 (70.0%) 17093 (30.0%) 57051 (100.0%)
Chi.squared df p.value
168.1 2 0
names(df)

[1] “MUNPAR.n” “SEMANAS” “CESAREA.f”
[4] “PESON1” “SEXO1.f” “NUMHV”
[7] “EDADM.2” “EDADP.2” “ESTUDIOM.f”
[10] “ESTUDIOP.f” “PAISNACM.n” “PAISNACP.n”
[13] “fecha_parto” “fecha_parto_anterior”

3.2 Tablas con janitor

La verdad es que janitor es un paquete fantástico!!!

janitor has simple functions for examining and cleaning dirty data. It was built with beginning and intermediate R users in mind and is optimized for user-friendliness. Advanced R users can already do everything covered here, but with janitor they can do it faster and save their thinking for the fun stuff.

  • hacer (y dar formato) a tablas de 1, 2 o 3 variables:
df %>% janitor::tabyl(CESAREA.f) %>% gt::gt()
CESAREA.f n percent
Sin cesárea 39958 0.7003909
Con cesárea 17093 0.2996091
df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f) %>% gt::gt()
CESAREA.f Primarios Medios Universidad
Sin cesárea 15734 11764 12460
Con cesárea 5769 5335 5989
df %>% janitor::tabyl(CESAREA.f, ESTUDIOM.f, SEXO1.f) #%>% gt::gt()
#> $Varón
#>    CESAREA.f Primarios Medios Universidad
#>  Sin cesárea      7969   6004        6499
#>  Con cesárea      3102   2752        3208
#> 
#> $Mujer
#>    CESAREA.f Primarios Medios Universidad
#>  Sin cesárea      7765   5760        5961
#>  Con cesárea      2667   2583        2781

3.3 Tablas con R-base

Las tablas con R-base tampoco están mal. El problema que tienen es que los resultados no se almacenan en data.frames:

table(df$CESAREA.f, df$ESTUDIOM.f, useNA = "always") 
#>              
#>               Primarios Medios Universidad  <NA>
#>   Sin cesárea     15734  11764       12460     0
#>   Con cesárea      5769   5335        5989     0
#>   <NA>                0      0           0     0
my_tabla <- table(df$CESAREA.f, df$ESTUDIOM.f) 
my_tabla %>% as.data.frame() %>% gt::gt()
Var1 Var2 Freq
Sin cesárea Primarios 15734
Con cesárea Primarios 5769
Sin cesárea Medios 11764
Con cesárea Medios 5335
Sin cesárea Universidad 12460
Con cesárea Universidad 5989
zz <- prop.table(my_tabla,1)
mosaicplot(t(zz), color = TRUE, main = "% de Cesáreas para niveles educativos de la madre")

3. Contrastes

Evidentemente R permite implementar múltiples técnicas y modelos estadísticos, así como hacer múltiples y variados contrastes de hipótesis. En esta pagina web tenéis un buena introducción a estos tema. Veamos algún ejemplo: por ejemplo:

t.test(df$PESON1)
t.test(df$PESON1 ~ df$CESAREA.f)
t.test(df$PESON1 ~ df$SEXO1.f) %>% psycho::analyze()
t.test(df$PESON1 ~ df$SEXO1.f, var.equal = TRUE, conf.level = .95) %>% psycho::analyze() %>% summary()
  • correlación entre dos variables cuantitativas
cor_results <- cor.test(df$PESON1, df$EDADM.2)  # Compute a correlation and store its result
psycho::analyze(cor_results)  # Run the analyze function on the correlation
#> The Pearson's product-moment correlation between df$PESON1 and df$EDADM.2 is not significant, very small, and positive (r(57049.00) = 0.00, 95% CI [-0.01, 0.01], p > .1).
  • un ejemplo para ver si hay diferencias en el peso en función de los estudios de la madre (!!!!)
library(purrr)
library(broom)
df %>% group_by(ESTUDIOM.f) %>% 
  summarise(t_test = list(t.test(PESON1))) %>% 
  mutate(tidied = map(t_test, tidy)) %>% 
  tidyr::unnest(tidied) %>% 
  ggplot(aes(estimate, ESTUDIOM.f)) + geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) + 
  labs(y = "")

chisq.test(df$CESAREA.f, df$SEXO1.f)
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  df$CESAREA.f and df$SEXO1.f
#> X-squared = 15.155, df = 1, p-value = 0.00009901
chisq.test(df$CESAREA.f, df$ESTUDIOM.f) 
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  df$CESAREA.f and df$ESTUDIOM.f
#> X-squared = 168.07, df = 2, p-value < 0.00000000000000022
  • Aquí puedes encontrar un curso sobre como implementar los contrastes estadísticos más habituales con R.

En este post, Jonas Kristoffer Lindeløv, nos presenta un cuadro con los contrastes de hipótesis más frecuentes y cómo efectuar esos contrastes en el contexto de los modelos de regresión con R.



4. Modelos

En este apartado vamos a ver como estimar modelos lineales y no lineales con R. Para una introducción a la estimación de modelos en R puedes ir aquí.

Tenemos/tengo suerte de que me puedo apoyar en lo que ya habéis visto en Econometría y otras asignaturas. El objetivo, aparte de ver como se pueden hacer análisis de regresión con R, es ir preparando el camino para introducir algo de Machine Learning.

Antes restrinjamos un poco el df:

df_m <- df %>% select(PESON1, SEMANAS, SEXO1.f, CESAREA.f, EDADM.2, EDADP.2, ESTUDIOM.f, ESTUDIOP.f)

4.1 Modelos lineales

La función para estimar modelos lineales es lm(), así que vamos a utilizarla para estimar nuestro primer modelo con R. Estimaremos un modelo lineal con variable a explicar el peso del bebe (PESON1) en función de todas las demás variables en df_m.

mod_1 <- lm(PESON1 ~ . , data = df_m)
summary(mod_1)

Generalmente las variables categóricas se introducen en los modelos mediante variables dummies. Crear dummies en R es sencillo, solo tienes que tener los datos como factor o como texto y R creará las dummies por ti cuando introduzcas la variable en lm(). Eso sí, es más facil elegir la categoría de referencia si la variable es un factor.

Para entender como fija los regresores para las variables categóricas, mira esto:

levels(df$SEXO1.f)
#> [1] "Varón" "Mujer"
levels(df$ESTUDIOM.f)
#> [1] "Primarios"   "Medios"      "Universidad"

Si necesitas cambiar la categoría de referencia siempre puedes usar forcast::fct_relevel()

zz <- forcats::fct_relevel(df$SEXO1.f, "Mujer")
levels(zz)
#> [1] "Mujer" "Varón"

Veamos que hay en el objeto mod_1

str(mod_1)
#listviewer::jsonedit(mod_2, mode = "view") ## Interactive option

Especificación del modelo

Lógicamente, a veces querremos seleccionar las variables explicativas:

mod_1 <- lm(PESON1 ~ SEMANAS, data = df)
mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + CESAREA.f + EDADM.2 , data = df_m)
mod_1 <- lm(log(PESON1) ~ log(SEMANAS), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = log(PESON1) ~ log(SEMANAS), data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.80542 -0.08174  0.00610  0.08778  1.39651 
#> 
#> Coefficients:
#>              Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)   0.63860    0.04359   14.65 <0.0000000000000002 ***
#> log(SEMANAS)  2.02921    0.01190  170.47 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1476 on 57049 degrees of freedom
#> Multiple R-squared:  0.3375, Adjusted R-squared:  0.3375 
#> F-statistic: 2.906e+04 on 1 and 57049 DF,  p-value: < 0.00000000000000022

Si queremos introducir interacciones entre los regresores, podemos usar el operador :, aunque casi mejor hacerlo directamente con I():

mod_1 <- lm(PESON1 ~ SEMANAS + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*EDADM.2), data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * EDADM.2), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2633.6  -279.6   -14.5   262.9  3709.3 
#> 
#> Coefficients:
#>                          Estimate   Std. Error t value             Pr(>|t|)    
#> (Intercept)          -2252.786335    37.656076  -59.83 < 0.0000000000000002 ***
#> SEMANAS                138.961398     0.994026  139.80 < 0.0000000000000002 ***
#> I(SEMANAS * EDADM.2)     0.060934     0.008869    6.87     0.00000000000647 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.2 on 57048 degrees of freedom
#> Multiple R-squared:  0.272,  Adjusted R-squared:  0.2719 
#> F-statistic: 1.066e+04 on 2 and 57048 DF,  p-value: < 0.00000000000000022

Si queremos introducir las variables originales y también las interacciones entre ellas, podemos hacerlo directamente o o utilizar el operador *:

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + SEMANAS:EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*EDADM.2, data = df_m)
mod_1 <- lm(PESON1 ~ SEMANAS*ESTUDIOM.f, data = df_m)

summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS * ESTUDIOM.f, data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2630.8  -280.8   -15.4   263.9  3727.1 
#> 
#> Coefficients:
#>                                 Estimate Std. Error t value            Pr(>|t|)
#> (Intercept)                   -2287.0303    57.9094 -39.493 <0.0000000000000002
#> SEMANAS                         141.5391     1.4842  95.366 <0.0000000000000002
#> ESTUDIOM.fMedios                -29.0793    90.4969  -0.321              0.7480
#> ESTUDIOM.fUniversidad           183.7079    91.1839   2.015              0.0439
#> SEMANAS:ESTUDIOM.fMedios          0.9445     2.3193   0.407              0.6838
#> SEMANAS:ESTUDIOM.fUniversidad    -4.0701     2.3383  -1.741              0.0818
#>                                  
#> (Intercept)                   ***
#> SEMANAS                       ***
#> ESTUDIOM.fMedios                 
#> ESTUDIOM.fUniversidad         *  
#> SEMANAS:ESTUDIOM.fMedios         
#> SEMANAS:ESTUDIOM.fUniversidad .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.2 on 57045 degrees of freedom
#> Multiple R-squared:  0.2719, Adjusted R-squared:  0.2718 
#> F-statistic:  4260 on 5 and 57045 DF,  p-value: < 0.00000000000000022

Recuerda que si queremos introducir algún regresor que sea la multiplicación de dos variables, tendremos que hacerlo con I()

mod_1 <- lm(PESON1 ~ SEMANAS + I(SEMANAS*SEMANAS), data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ SEMANAS + I(SEMANAS * SEMANAS), data = df_m)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -2635.6  -282.7   -15.6   261.5  3803.8 
#> 
#> Coefficients:
#>                        Estimate Std. Error t value             Pr(>|t|)    
#> (Intercept)          -3116.4624   255.2162 -12.211 < 0.0000000000000002 ***
#> SEMANAS                188.5366    13.8452  13.617 < 0.0000000000000002 ***
#> I(SEMANAS * SEMANAS)    -0.6514     0.1878  -3.469             0.000522 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 434.3 on 57048 degrees of freedom
#> Multiple R-squared:  0.2715, Adjusted R-squared:  0.2715 
#> F-statistic: 1.063e+04 on 2 and 57048 DF,  p-value: < 0.00000000000000022

Cambien puede sernos de utilidad la función poly()

lm(PESON1 ~ poly(SEMANAS, degree=3), data = df_m)
#> 
#> Call:
#> lm(formula = PESON1 ~ poly(SEMANAS, degree = 3), data = df_m)
#> 
#> Coefficients:
#>                (Intercept)  poly(SEMANAS, degree = 3)1  
#>                       3238                       63314  
#> poly(SEMANAS, degree = 3)2  poly(SEMANAS, degree = 3)3  
#>                      -1507                       -7131

Resultados de estimación

Una vez que sabemos la sintaxis de stats::lm(), vamos a ver como podemos acceder a la información que devuelve lm(). ya hemos visto que los resultados se almacenan en una lista, así que podemos utilizar los operadores habituales de las listas para acceder a la información. Por ejemplo:

zz_betas     <- mod_1[[1]]

zz_residuals <- mod_1[[2]]
zz_residuals <- mod_1[["residuals"]]
zz_residuals <- mod_1$residuals

zz_betas_2 <- mod_1[1]                #- veis la diferencia
zz_betas_3 <- mod_1[[1]]              #- veis la diferencia

zz_betas_2a <- mod_1["coefficients"]
zz_betas_3a <- mod_1$"coefficients"

Afortunadamente ya hay construidas algunas funciones para manipular/ver los resultados de estimación. Por ejemplo:

library(ggfortify)
summary(mod_1)                   #- tabla resumen
summary(mod_1)$coefficients      #- tabla resumen con los coeficientes
coefficients(mod_1)              #- coeficientes estimados
confint(mod_1, level = 0.95)     #- Intervalos de confianza para los coeficientes
#fitted(mod_1)                   #- predicciones (y-sombrero, y-hat)
#residuals(mod_1)                #- vector de residuos
#model.matrix(mod_1)             #- extract the model matrix
anova(mod_1)                     #- ANOVA
vcov(mod_1)                      #- matriz de var-cov para los coeficientes 
#influence(mod_1)                #- regression diagnostics 
# diagnostic plots
layout(matrix(c(1,2,3,4),2,2))   #- optional 4 graphs/page
plot(mod_1)                      #-
library(ggfortify)
autoplot(mod_1, which = 1:6, ncol = 2, colour = "steelblue")

Para hacer predicciones con observaciones fuera de la muestra has de usar la función predict, proporcionándole las observaciones a predecir como un df.

nuevas_observaciones <- df_m %>% slice(c(3,44,444))
predict(mod_1, newdata = nuevas_observaciones)  #- predicciones puntuales
predict(mod_1, newdata = nuevas_observaciones, type = 'response', se.fit = TRUE)  #- tb errores estándar  predictions
predict(mod_1, newdata = nuevas_observaciones, interval = "confidence")  #- intervalo (para el valor esperado)
predict(mod_1, newdata = nuevas_observaciones, interval = "prediction")  #- intervalo (para valores individuales)

Hay muchas más funciones para valorar la idoneidad de un modelo. Por ejemplo aquí.

El paquete GGally permite hacer:

  • gráfico de los coeficientes estimados
GGally::ggcoef(mod_1)

df_jugar <- df_m %>% select(EDADM.2, SEXO1.f, SEMANAS)
GGally::ggpairs(df_jugar)

Errores estándar robustos

Lo habitual es que al estimar un modelo, tengamos situaciones de heterocedasticidad, clustering etc … Podemos ajustar los errores estándar a estas situaciones. El paquete de referencia solía ser sandwich, aunque quizás ahora el paquete de referencia sea estimatr. Por ejemplo se pueden obtener errores estándar robustos con estimatr::lm_robust()2.

#- install.packages("emmeans")
mod_1 <- lm(PESON1 ~ SEMANAS, data = df_m)
mod_1_ee <- estimatr::lm_robust(PESON1 ~ SEMANAS, data = df_m)
broom::tidy(mod_1_ee, mod_1, conf.int = T)
#>          term   estimate std.error statistic
#> 1 (Intercept) -2240.6797 60.139009 -37.25834
#> 2     SEMANAS   140.6182  1.533625  91.69008
#>                                                                                                                                                                                                                                                                                                                p.value
#> 1 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003194318
#> 2 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
#>     conf.low  conf.high    df outcome
#> 1 -2358.5525 -2122.8069 57049  PESON1
#> 2   137.6123   143.6242 57049  PESON1

Paquete broom

Un opción interesante es utilizar el paquete broom. Este paquete tiene 3 funciones útiles:

  • tidy()

  • augment()

  • glance()

zz <- broom::tidy(mod_1, conf.int = TRUE)
zz %>% gt::gt()
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -2240.6797 37.6300489 -59.54496 0 -2314.4348 -2166.9246
SEMANAS 140.6182 0.9647186 145.76088 0 138.7274 142.5091
mod_1 %>% broom::glance() %>% select(adj.r.squared, p.value)
#> # A tibble: 1 x 2
#>   adj.r.squared p.value
#>           <dbl>   <dbl>
#> 1         0.271       0
broom::glance(mod_1)
#> # A tibble: 1 x 11
#>   r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC    BIC
#>       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>  <dbl>
#> 1     0.271         0.271  434.    21246.       0     2 -4.27e5 8.55e5 8.55e5
#> # ... with 2 more variables: deviance <dbl>, df.residual <int>
zz <- broom::augment(mod_1)
  • un ejemplo de utilidad de broom
mod_1 %>% tidy() %>% filter(p.value < 0.05)
#> # A tibble: 2 x 5
#>   term        estimate std.error statistic p.value
#>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 (Intercept)   -2241.    37.6       -59.5       0
#> 2 SEMANAS         141.     0.965     146.        0
  • Un ejemplo de uso de broom, pero antes vamos a recordar alguna cosa de ggplot2:
ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = ESTUDIOM.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
      geom_point(alpha = 0.1) +  geom_smooth(method = "lm")

ggplot(data = df_m, mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
      geom_point(alpha = 0.1) +  geom_smooth()

Ahora sí viene el ejemplo en que se usa el paquete broom

mod_1 <- lm(PESON1 ~ EDADM.2 + SEXO1.f , data = df_m)
summary(mod_1)
#> 
#> Call:
#> lm(formula = PESON1 ~ EDADM.2 + SEXO1.f, data = df_m)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -2841.04  -280.59    19.01   318.31  3058.89 
#> 
#> Coefficients:
#>                Estimate Std. Error t value            Pr(>|t|)    
#> (Intercept)  3288.96818   13.29270 247.427 <0.0000000000000002 ***
#> EDADM.2         0.07408    0.40164   0.184               0.854    
#> SEXO1.fMujer -110.83978    4.23833 -26.152 <0.0000000000000002 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 505.8 on 57048 degrees of freedom
#> Multiple R-squared:  0.01185,    Adjusted R-squared:  0.01181 
#> F-statistic:   342 on 2 and 57048 DF,  p-value: < 0.00000000000000022

td_mod_1 <- mod_1 %>% broom::augment(data = df_m) 

td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = PESON1,  color = SEXO1.f)) +
                geom_point(alpha = 0.1) +
                geom_line(aes(y = .fitted, group = SEXO1.f))

td_mod_1 %>% ggplot(mapping = aes(x = EDADM.2, y = .fitted,  color = SEXO1.f)) +
                geom_line(aes(group = SEXO1.f))

(!!!!) Por ejemplo también permite fácilmente estimar modelos por grupos:

df_m %>% group_by(SEXO1.f) %>% do(tidy(lm(PESON1 ~ SEMANAS, .)))
#> # A tibble: 4 x 6
#> # Groups:   SEXO1.f [2]
#>   SEXO1.f term        estimate std.error statistic p.value
#>   <fct>   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#> 1 Varón   (Intercept)   -2363.     52.1      -45.4       0
#> 2 Varón   SEMANAS         145.      1.34     109.        0
#> 3 Mujer   (Intercept)   -2158.     53.2      -40.5       0
#> 4 Mujer   SEMANAS         137.      1.36     100.        0

O hacer gráficos de los intervalos de los coeficientes:

td_mod_1 <- tidy(mod_1, conf.int = TRUE)
ggplot(td_mod_1, aes(estimate, term, color = term)) +
    geom_point() +
    geom_errorbarh(aes(xmin = conf.low, xmax = conf.high)) +
    geom_vline(xintercept = 0)

o este:

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line()

mod_1 %>% augment(data = df_m) %>%
 ggplot(mapping = aes(x = SEMANAS , y = .fitted, color = SEXO1.f)) +
   geom_point(mapping = aes(y = PESON1), alpha = 0.1) +
   geom_line() +
  facet_wrap(vars(ESTUDIOM.f))


Comparación de modelos

mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2,  data = df_m)
mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)

anova(mod_1, mod_2)
#> Analysis of Variance Table
#> 
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#>   Res.Df         RSS Df Sum of Sq      F  Pr(>F)  
#> 1  57048 10755232331                              
#> 2  57047 10754104306  1   1128025 5.9838 0.01444 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(mod_1, mod_2)
#>       df      AIC
#> mod_1  4 854907.8
#> mod_2  5 854903.8
lmtest::lrtest(mod_1, mod_2)    
#> Likelihood ratio test
#> 
#> Model 1: PESON1 ~ SEMANAS + EDADM.2
#> Model 2: PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS * EDADM.2)
#>   #Df  LogLik Df  Chisq Pr(>Chisq)  
#> 1   4 -427450                       
#> 2   5 -427447  1 5.9839    0.01444 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(!!!) Vamos a ordenar los modelos en función de su AIC. Para ello vamos a crear una lista con modelos

modelos <- list(mod_1 <- lm(PESON1 ~ SEMANAS + EDADM.2,  data = df_m),
                mod_2 <- lm(PESON1 ~ SEMANAS + EDADM.2 + I(SEMANAS*EDADM.2), data = df_m)  )

modelos_ordered_AIC <- purrr::map_df(modelos, broom::glance, .id = "model") %>% arrange(AIC)

modelos_ordered_AIC %>% gt::gt()
model r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual
2 0.2720188 0.2719805 434.1809 7105.43 0 4 -427446.9 854903.8 854948.6 10754104306 57047
1 0.2719425 0.2719169 434.1999 10654.22 0 3 -427449.9 854907.8 854943.6 10755232331 57048
  • una tabla con broom (antes creamos una función (!!!!)):
my_kable <- function(df){ gt::gt(mutate_if(df, is.numeric, round, 2)) }

tidy(mod_1) %>% my_kable
term estimate std.error statistic p.value
(Intercept) -2327.62 39.76 -58.55 0
SEMANAS 140.92 0.97 145.97 0
EDADM.2 2.33 0.35 6.75 0

Con el paquete modelr podemos fácilmente comparar las predicciones de varios modelos:

zz <- df_m %>% modelr::gather_predictions(mod_1, mod_2)

Ejercicio: Pasad zz a formato ancho para ver las predicciones de los 2 modelos:

zz1 <- pivot_wider(zz, names_from = model, values_from = pred)

4.2 Modelos GLM

Por ejemplo un Logit:

mod_logit <- glm(CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f, family = binomial(link = "logit"), data = df_m)

summary(mod_logit)
#> 
#> Call:
#> glm(formula = CESAREA.f ~ PESON1 + SEMANAS + EDADM.2 + ESTUDIOM.f, 
#>     family = binomial(link = "logit"), data = df_m)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -1.8803  -0.8662  -0.7599   1.3804   2.1197  
#> 
#> Coefficients:
#>                          Estimate  Std. Error z value             Pr(>|z|)    
#> (Intercept)            2.23616479  0.20327735  11.001 < 0.0000000000000002 ***
#> PESON1                 0.00013040  0.00002128   6.129       0.000000000882 ***
#> SEMANAS               -0.13285957  0.00567450 -23.413 < 0.0000000000000002 ***
#> EDADM.2                0.04918850  0.00191326  25.709 < 0.0000000000000002 ***
#> ESTUDIOM.fMedios       0.10281855  0.02322645   4.427       0.000009564693 ***
#> ESTUDIOM.fUniversidad  0.07570366  0.02335789   3.241              0.00119 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 69663  on 57050  degrees of freedom
#> Residual deviance: 68142  on 57045  degrees of freedom
#> AIC: 68154
#> 
#> Number of Fisher Scoring iterations: 4

En los modelos lineales, calcular efectos marginales es sencillo, pero el Logit es un modelo no lineal. Para calcular efectos marginales en modelos no lineales podemos usar el paquete margins. Por ejemplo, calculemos el Efecto marginal medio o average marginal effect (AME) en mod_logit.

mod_1_AME <- mod_logit %>% margins::margins() %>% summary() 
mod_1_AME
#>                 factor     AME     SE        z      p   lower   upper
#>                EDADM.2  0.0100 0.0004  26.1951 0.0000  0.0093  0.0108
#>       ESTUDIOM.fMedios  0.0210 0.0047   4.4201 0.0000  0.0117  0.0303
#>  ESTUDIOM.fUniversidad  0.0154 0.0047   3.2374 0.0012  0.0061  0.0247
#>                 PESON1  0.0000 0.0000   6.1353 0.0000  0.0000  0.0000
#>                SEMANAS -0.0271 0.0011 -23.8118 0.0000 -0.0294 -0.0249

Si queremos calcular efectos marginales para unos valores concretos de los regresores

mod_logit %>% margins::margins(at = list(SEMANAS = c(25, 35), ESTUDIOM.f = c("Primarios", "Medios", "Universidad")),  variables = "EDADM.2" ) 

Si prefieres visualizarlo, utiliza margins::cplot():

margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2", what = "effect", drop = TRUE)
margins::cplot(mod_logit, x = "ESTUDIOM.f", dx = "EDADM.2")


4.3 Tablas para modelos

Con sjPLot y friends …

m1 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2, data = df)
m2 <- lm(PESON1 ~ SEMANAS + SEXO1.f + EDADM.2 + EDADP.2 + ESTUDIOM.f + ESTUDIOP.f, data = df)
sjPlot::tab_model(m1)
sjPlot::plot_model(m1, sort.est = TRUE)
sjPlot::tab_model(m1, m2)


Con stargazer

mod_1 <- lm(PESON1 ~ SEMANAS + SEXO1.f , data = df_m)
stargazer::stargazer(mod_1, type = "html")
Dependent variable:
PESON1
SEMANAS 141.401***
(0.955)
SEXO1.fMujer -123.584***
(3.604)
Constant -2,211.576***
(37.258)
Observations 57,051
R2 0.286
Adjusted R2 0.286
Residual Std. Error 429.964 (df = 57048)
F Statistic 11,430.040*** (df = 2; 57048)
Note: p<0.1; p<0.05; p<0.01


Con modelsummary

#remotes::install_github('vincentarelbundock/modelsummary')
library(modelsummary)

mys_modelitos <- list()
mys_modelitos[["PESON1:  OLS 1"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1.f,            df_m)
mys_modelitos[["PESON1:  OLS 2"]]    <-  lm(PESON1    ~ SEMANAS + SEXO1.f + EDADM.2 , df_m)
mys_modelitos[["CESAREA: Logit 1"]]  <- glm(CESAREA.f ~ SEMANAS + SEXO1.f , data = df_m, family = binomial(link = "logit"))

mm <- msummary(mys_modelitos, title = "Resultados de estimación")
mm
Resultados de estimación
PESON1: OLS 1 PESON1: OLS 2 CESAREA: Logit 1
(Intercept) -2211.576 -2301.307 3.761
(37.258) (39.360) (0.184)
SEMANAS 141.401 141.711 -0.118
(0.955) (0.956) (0.005)
SEXO1.fMujer -123.584 -123.749 -0.061
(3.604) (3.602) (0.018)
EDADM.2 2.406
(0.342)
Num.Obs. 57051 57051 57051
R2 0.286 0.287
Adj.R2 0.286 0.287
AIC 853789.1 853741.5 69027.1
BIC 853824.9 853786.2 69054.0
Log.Lik. -426890.542 -426865.740 -34510.573


Con reports

  • informes con el paquete reports
library(reports) #- devtools::install_github("neuropsychology/report")
my_model <- lm(PESON1 ~ SEMANAS + SEXO1.f, df_m)
rr <- report(my_model, target = 1)
rr


to_table(rr) %>% gt::gt()


to_fulltext(rr)


Recuerda también que se puede obtener la ecuación de un modelo con:

library(equatiomatic)  #- remotes::install_github("datalorax/equatiomatic")
extract_eq(mod_1)
#> $$
#> \text{PESON1} = \alpha + \beta_{1}(\text{SEMANAS}) + \beta_{2}(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$
extract_eq(mod_1, use_coefs = TRUE)
#> $$
#> \text{PESON1} = -2211.58 + 141.4(\text{SEMANAS}) - 123.58(\text{SEXO1.f}_{\text{Mujer}}) + \epsilon
#> $$



5. Otros modelos/técnicas

Para que los resultados salgan rápido y se puedan visualizar, vamos a restringir el df:

df_x <- df %>% select(SEMANAS, CESAREA.f, EDADM.2, PESON1, SEXO1.f) %>%  #- cogemos unas pocas variables 
                na.omit() %>%                                            #- quitamos NA's
                dplyr::sample_n(200) %>%                                 #- y solo 200 bebes
                mutate_if(is.factor, as.numeric)                         #- tienen q ser numericas

Componentes principales

  • usamos la función DataExplorer::plot_prcomp() para visualizar los resultados obtenidos con stats::prcomp()
DataExplorer::plot_prcomp(df_x, variance_cap = 0.9, nrow = 2L, ncol = 2L)
  • FactoMineR::PCA()
FactoMineR::PCA(df_x)
  • El paquete factoextra para visualizar
res.pca <- FactoMineR::PCA(df_x,  ncp = 5, graph = FALSE)
factoextra::fviz_screeplot(res.pca, addlabels = TRUE)
factoextra::fviz_pca_var(res.pca, col.var = "black")
factoextra::fviz_contrib(res.pca, choice = "var", axes = 1, top = 10) # Contributions of variables to PC1
factoextra::fviz_pca_ind(res.pca)   #- resultados para bebes individuales
factoextra::fviz_pca_ind(res.pca, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)# Avoid text overlapping (slow if many points))

  • psycho::n_factors(): Find optimal components number using maximum method aggreement, pero no me acaba de funcionar
library(psycho)
results <- df_numeric %>% psycho::n_factors()
print(results)
plot(results)
summary(correlation(df_numeric))


Cluster

df_xx <- base::scale(df_x)
  • Cluster:
k3 <- stats::kmeans(df_x, centers = 3, nstart = 7) 
factoextra::fviz_cluster(k3, data = df_x)

Hice este gráfico para ver si lo entendía, pero …

df_x %>% as_tibble() %>% mutate(cluster = k3$cluster) %>%
  ggplot(aes(SEMANAS, PESON1, shape = as.factor(CESAREA.f), color = factor(cluster), label = SEXO1.f)) +
  geom_point()  +
  geom_text()
  • Varios clusters:
#- https://uc-r.github.io/kmeans_clustering
k2 <- kmeans(df_x, centers = 2, nstart = 7)
k3 <- kmeans(df_x, centers = 3, nstart = 7)
k4 <- kmeans(df_x, centers = 4, nstart = 7)
k5 <- kmeans(df_x, centers = 5, nstart = 7)
p1 <- fviz_cluster(k2, geom = "point",  data = df_x) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = df_x) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = df_x) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = df_x) + ggtitle("k = 5")
library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)

#- https://uc-r.github.io/kmeans_clustering
distance <- factoextra::get_dist(df_x)
factoextra::fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
  • número óptimo de clusters:
factoextra::fviz_nbclust(df_xx, kmeans, method = "gap_stat")


Arboles de decisión

  • Arboles de decisión con la función explore::explain_tree() que: Explain a target using a simple decision tree (classification or regression)
df_m %>% explore::explain_tree(target = PESON1)
df_m %>% explore::explain_tree(target = CESAREA.f)
df_m %>% explore::explain_tree(target = SEXO1.f)
df_m %>% explore::explain_tree(target = ESTUDIOM.f)
fit <- rpart::rpart(SEXO1.f ~ . , data = df_m)
rattle::fancyRpartPlot(fit, palettes = c("Greens", "Reds"), sub = "")




6. ML

Simplemente algunas referencias:


  1. El warning nos salta porque algunas variables son factores y tienen NA’s, pero NA no es uno de los levels del factor. El warning simplemente nos recuerda que hay como un level “oculto”.

  2. Por defecto, el paquete usa Eicker-Huber-White robust standard errors, habitualmente conocidos como errores estándar “HC2”. Se pueden especificar otros métodos con la opción se_type. Por ejemplo se puede utilizar el método usado por defecto en Stata. Aquí puedes encontrar porque los errores estándar de Stata difieren de los usados en R y Phyton